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Abstract 

This  paper  introduces  a  new  approach  to  the  problem  of  quantitatively  reconstructing 
cylindrically  symmetric  objects  from  radiograph  data  obtained  via  x-ray  tomography.  Specif¬ 
ically,  a  mixed  variable  programming  (MVP)  problem  is  formulated,  in  which  the  variables 
of  interest  are  the  number  and  types  of  materials  and  the  thickness  of  each  concentric  layer. 

The  objective  function  is  a  measure  of  distance  between  one-dimensional  radiograph  data 
and  a  material  property  vector  operated  on  by  a  forward  projection  based  on  the  Abel  trans¬ 
form.  The  mixed  variable  pattern  search  (MVPS)  algorithm  for  linearly  constrained  MVP 
problems  is  applied  to  the  problem  by  means  of  the  NOMADm  MATLAB  software  package. 
Numerical  results  are  presented  for  several  test  configurations  and  show  that,  while  there 
are  difficulties  yet  to  be  overcome,  the  method  appears  to  be  very  promising  for  solving  this 
class  of  problems  in  practice. 
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1  Introduction 

Tomography  refers  to  the  “cross-sectional  imaging  of  an  object  from  either  transmission  or 
reflection  of  data  collected  by  illuminating  the  object  from  many  different  directions”  [211 .  While 
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the  majority  of  research  in  computerized  tomographic  techniques  has  been  focused  on  diagnostic 
medicine  ( e.g .,  radioisotopes,  ultrasound,  magnetic  resonance),  there  has  also  been  significant 
research  in  other  areas,  such  as  mapping  of  underground  resources,  nondestructive  testing  of 
engineered  parts  such  as  rocket  engines,  brightness  distribution  determination  of  a  celestial 
sphere,  and  three-dimensional  imaging  using  an  electron  microscope  pi]. 

Mathematical  inversion  of  cylindrically  symmetric  objects  was  originally  solved  analytically 
by  Abel  [2]  in  1826,  resulting  in  Abel  inversion  techniques.  In  1917,  Radon  discovered  a  way  to 
mathematically  reconstruct  any  function  from  its  projections,  but  his  methods  produce  inexact 
reconstructions  that  are  only  useful  in  qualitative  analysis,  such  as  medical  imaging.  In  1972 
Houndsfield  invented  the  first  x-ray  computed  tomographic  scanner,  which  could  reconstruct  an 
object’s  image  from  its  x-ray  projections  |21j. 

In  this  paper,  we  consider  the  problem  of  quickly  and  quantitatively  reconstructing  objects 
from  x-ray  radiograph  data.  In  particular,  we  make  use  of  the  Abel  transform  p3j  and  mixed 
variable  optimization  mm  to  determine  the  composition  of  a  cylindrically  symmetrical  object 
made  of  concentric  material  layers,  in  which  the  composition  is  defined  by  the  number,  types, 
and  thicknesses  of  its  layers.  This  class  of  problems  is  challenging  because  the  material  types 
are  nonnumeric  categorical  variables. 

Figure  [l]  shows  the  lateral  section  of  a  cylindrically  symmetric  object  composed  of  four 
material  layers.  The  shading  represents  different  materials,  while  the  concentric  circles  denote 
the  material  interfaces.  The  dashed  lines  represent  the  x-rays  as  they  travel  through  the  object, 
while  the  horizontal  line  at  the  bottom  of  the  picture  represents  the  detector  used  to  record 
radiograph  data. 

Current  approaches  involve  either  the  inversion  of  the  Abel  transform  or  some  type  of  regu¬ 
larization  scheme.  However,  both  approaches  have  inherent  difficulties  and  have  not  performed 
well  for  quantitative  reconstructions.  As  an  alternative,  we  introduce  a  new  approach  in  which 
the  problem  is  formulated  as  a  mixed  variable  optimization  problem  and  solved  numerically 
using  the  Audet-Dennis  mixed  variable  generalized  pattern  search  (MVPS)  algorithm  [11J.  In 
this  case,  the  objective  is  to  minimize  a  measure  of  error  between  x-ray  radiograph  data  and 
possible  object  configurations.  Mixed  variable  problems  are  those  that  can  have  both  continuous 
and  categorical  variables,  the  latter  being  discrete  variables  that  must  always  take  their  values 
from  a  predefined  list  or  set. 

The  remainder  of  the  paper  is  laid  out  as  follows.  In  Section  [2]  we  describe  the  x-ray  to¬ 
mography  problem  in  detail  and  give  the  mixed  variable  formulation  of  the  problem.  Section  [3] 
describes  the  generalized  pattern  search  algorithm  for  linearly  constrained  mixed  variable  opti¬ 
mization  problems.  Section  [4]  describes  the  experimental  conditions  for  numerical  testing,  and 
Section  [5]  presents  computational  results  for  several  test  configurations.  Finally,  Section  [6] offers 
some  concluding  remarks. 

2  The  Quantitative  Object  Reconstruction  Problem 

X-rays  are  electromagnetic  waves  propagating  through  space,  forming  beams  of  photons.  When 
these  beams  interact  with  electrons  or  protons  within  the  atoms  of  a  material,  they  are  affected 
by  attenuation  (loss  of  energy)  and  photon  scattering. 
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Figure  1:  Cylindrically  Spherical  Object  being  x-rayed. 


At  the  simplest  level,  an  x-ray  radiograph  is  a  projection  of  x-ray  beam  intensity  as  the 
beams  pass  through  a  certain  amount  of  material  between  the  x-ray  source  and  the  detector.  If 
several  radiographs  are  taken  from  different  viewing  angles,  it  may  be  possible  to  reconstruct 
the  shape  of  the  original  object. 

In  a  linear  attenuation  model,  the  intensity  I  of  the  x-ray  beam  is  related  to  the  attenuation 
properties  of  the  material  by  the  equation, 

I  =  h  exp  J  n{r)d(^j  ,  (1) 

where  Iq  is  the  initial  intensity  of  the  x-ray  beam  and  //(r)  is  the  attenuation  coefficient,  which 
is  a  function  of  the  radial  distance  r  from  the  center  of  the  object.  The  line  integral  is  computed 
as  shown  in  Figure  [lj 

Reconstruction  of  objects  usually  requires  the  use  of  a  Radon  transform,  which  reduces 
to  the  Abel  transform  in  the  case  of  a  single  x-ray  radiograph  of  a  cylindrically  symmetric 
object.  However,  difficulties  arise  from  its  use  on  this  problem,  with  respect  to  both  the  physics 
of  radiography  and  the  properties  of  the  Abel  transform.  Not  only  do  multiple  material  layers 
complicate  attenuation,  but  there  are  also  inherent  difficulties  caused  by  photon  scattering  mm- 
The  combination  of  all  of  these  processes  makes  quantitative  reconstruction  difficult  |8i . 

Another  difficulty  comes  from  the  type  of  x-ray  device  used.  In  many  applications,  a  poly¬ 
chromatic  x-ray  source  must  be  used.  However,  a  polychromatic  beam  exacerbates  the  problem 
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because  both  scattering  and  attenuation  are  energy  dependent  effects.  This  also  leads  to  a  phe¬ 
nomenon  known  as  beam  hardening,  in  which  the  beam’s  spectral  content  (the  distribution  of 
the  x-ray  beam’s  energy)  is  different  at  the  detector  than  at  the  source  [2T] . 

Despite  the  difficulties,  the  Abel  transform  remains  a  key  tool  in  describing  the  tomographic 
measurement  operator.  In  simplest  terms,  this  transform  is  a  linear  operator  $  that  maps  the 
object’s  properties  /r  to  the  measured  quantities  d  obtained  from  the  radiograph;  i.e. , 


<h/r  =  d, 


(2) 


which  ideally  could  be  solved  simply  by  taking  the  inverse  T-1.  However,  x-ray  detectors 
measure  photon  intensity,  and  in  order  to  determine  the  object’s  attenuation  properties,  intensity 
data  must  be  converted  to  an  areal  attenuation.  In  the  simple  linear  model  described  by  ([!]), 
d(x)  =  -log (7//0).  In  [El  E],  the  authors  consider  a  one-dinrensional  reconstruction  of  the 
object  from  a  single  radiograph  with  a  single  property  description,  /r(r).  Neglecting  the  physical 
interactions  of  scattering,  the  continuous  inverse  Abel  transform  for  this  case  may  be  expressed 
as 


1  d  f00  xd(x) 
vrr  dr  Jr  yjx2  -  r2  ’ 


(3) 


where  the  attenuation  coefficient  n  is  a  line  integral  relative  to  the  intensity  d  and  x  is  the 
distance  along  the  radiograph. 

However,  application  of  0  is  problematic  because  of  the  singularity  at  the  lower  limit  of 
integration.  Asaki  et  al.  |10|  suggest  discretizing  the  Abel  transform, 


d(x)  =  2 


rKr) 


dr, 


(4) 


and  then  solving  the  resulting  non-sparse  linear  system  P[i  =  d,  where  P  E  Mmxn,  /r  E  Mn,  and 
d  E  Mm.  This  procedure  is  a  common  practice  m-  While  discretization  avoids  the  integral 
singularity  in  it  may  still  produce  poor  results  for  several  reasons  mm-  First,  Q  is  a 
significantly  simplified  model  that  neither  accounts  for  the  physical  interactions  of  scattering 
nor  polychromatic  effects.  Second,  if  m  <  n,  P  is  not  invertible,  and  even  when  n  =  m  and  P 
exists,  P  is  ill-conditioned. 

To  overcome  these  difficulties,  Asaki  (7)  suggests  various  regularization  techniques  to  sup¬ 
press  noise  growth  caused  by  the  ill-conditioning  of  P.  Specifically,  he  proposes  solving  the 
optimization  problem, 


mm  F(fi)  =  || Pn  -  d||  +  aR(n),  (5) 

where  ||  ■  ||  is  a  data  fidelity  norm  that  measures  how  well  the  object  matches  the  data,  a  G  Mis  the 
regularization  parameter  chosen  to  estimate  the  variance  of  the  data,  and  R  is  the  regularization 
function.  Regularization  has  proven  to  be  successful  in  suppressing  noise,  but  it  is  unable 
to  incorporate  any  significant  prior  knowledge  about  the  object,  and  there  is  still  no  way  to 
quantitatively  determine  the  object’s  material  composition  except  in  the  simplest  cases  |9j. 

A  variety  of  other  methods  have  been  employed  to  overcome  the  difficulties  associated  with 
the  inverse  transform.  Basic  filtered  and  unfiltered  inversions  are  described  by  Dasch  m- 
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Polynomial  interpolation  methods  have  been  successful  in  describing  objects  of  smoothly  varying 
material  properties  pHUE].  Basis  function  expansions  are  another  popular  area;  see  [201  j!2?]  for 
some  examples.  Transform  techniques  are  demonstrated  by  Smith  et  al.  |29] . 

For  objects  of  piecewise  constant  material  properties,  the  most  interesting  solution  method 
is  that  of  Deutsch  et  al.  [18j  [19] .  They  provide  a  rapidly  converging  algorithm  for  determining 
the  attenuation  coefficients  of  a  layered  object.  The  method  requires  no  prior  knowledge  of 
attenuation  values,  but  does  assume  known  and  fixed  layer  boundary  positions. 

As  an  alternative,  we  formulate  the  problem  as  a  mixed  variable  optimization  problem  and 
apply  the  Audet-Dennis  mixed  variable  pattern  search  algorithm  to  solve  it.  Specifically,  we 
consider  objects  with  n  G  {nmjn, . . . ,  nmax}  concentric  material  layers,  each  with  a  material  type 
nii  £  M  and  outer  edge  location  Xi  >  0,  i  =  1, 2, . . . ,  n,  that  comes  from  a  library  (or  list)  of 
materials  M  and  has  a  minimum  thickness  of  6  >  0.  Given  a  radiograph  of  fixed  length  L  >  0 
(from  the  object  center  xo  to  the  edge  of  the  radiograph)  with  pixel  size  w  >  0,  radiograph  data 
d  of  size  nd,  and  discretized  Abel  Transform  matrix  P,  we  consider  the  optimization  problem, 


min 

n,m,x 


f(n,m,x) 


-^P/d(n,  m,  x)  —  ln(d) 
wz 


(6) 


s.  t.  5  <  Xi  —  Xi- 1  <  L  —  nd,  i  =  1, . . . ,  n, 
n5  <  xn  <  L  —  6, 


(7) 

(8) 


where  /r(n,  m,  x)  is  a  vector  of  material  properties  that  depend  on  the  number  of  material  layers 
n,  the  material  types  m  £  Mn,  and  material  edge  locations  x  £  Mn.  The  linear  constraints 
in  0-©  enforce  the  minimum  thickness  requirement  on  each  material  layer.  This  mixed 
variable  problem  is  particularly  challenging  and  interesting  because  the  categorical  variables  m 
are  nonnunreric  and  the  problem  dimension  n  is  itself  a  variable  in  the  problem. 

Our  formulation  as  a  mixed  variable  optimization  problem  differs  from  the  approach  of 
Deutsch,  et  al.  mm  in  several  ways.  First,  we  allow  both  the  number  and  position  of  material 
boundaries  to  vary.  This  allows  for  the  very  real  possibility  in  applications  that  these  quantities 
are  unknown  or  only  approximately  known.  Second,  our  proposed  method  does  not  rely  on 
a  linearized  transformation  between  the  object  description  and  the  radiograph  data.  This  is 
important  for  applications  in  which  the  imaging  particles  have  significant  nonlinear  interaction 
with  the  object  being  imaged.  Third,  our  method  provides  an  unabiguous  object  material 
description,  rather  than  a  scalar  attenuation  coefficient  (or  scalar  mass  density)  from  which  a 
material  must  be  inferred. 


3  Mixed  Variable  Pattern  Search 

Audet  and  Dennis  m  introduced  the  class  of  mixed  variable  pattern  search  (MVPS)  algorithms 
as  an  extension  of  the  original  pattern  search  algorithms  ESUaiEI]  and  proved  convergence 
to  appropriately  defined  [26]  stationary  points.  A  hierarchy  of  convergence  results,  depending 
on  the  smoothness  of  the  objective  function  is  given  in  j3).  Kokkolaras,  Audet  and  Dennis  [22] 
applied  the  MVPS  algorithm  to  the  design  of  a  thermal  insulation  system  and  showed  a  65% 
reduction  in  the  objective  function  value  over  previous  results  that  optimized  only  with  respect 
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to  the  continuous  variables.  MVPS  has  also  been  extended  to  problems  with  stochastic  noise 
in  the  objective  function  m  and  to  problems  with  nonlinear  constraints  ElEJ-  a  more  general 
framework  for  derivative-free  mixed  variable  optimization  is  described  in  |26j. 

The  class  of  MVPS  algorithms  is  designed  to  numerically  solve  optimization  problems  in 
which  each  variable  x  =  ( xc ,  xd)  is  partitioned  into  its  continuous  and  categorical  parts,  xc  E 
Xc  C  MnC  and  xd  E  Xd  C  Zn  ,  respectively,  where  nc  and  nd  denote  the  maximum  respective 
dimensions  of  these  variables.  We  adopt  the  convention  of  ignoring  unused  variables.  The 
general  MVP  problem  can  now  be  expressed  as 

min  f(x),  (9) 

where  /  :  X  — »  EU  {00} ,  and  the  domain  X  is  the  union  of  continuous  domains  across  possible 
categorical  variable  values;  i.e., 


x  =  U  (xV)  x  {A), 

xd£Xd 

with  the  convention  that  X  =  Xc  if  nd  =  0.  Furthermore,  Xc  is  defined  by  a  finite  set  of  bound 
and  linear  constraints,  dependent  on  the  values  of  xd.  That  is, 

Xc(xd )  =  {xc  E  Rn‘  :  £(xd )  <  A(xd)xc  <  u(xd )}, 

where  A(xd)  E  Wnxn'c  is  a  real  matrix,  £(xd),u(xd)  6  (RU  {±oo})m,  and  £(xd)  <  u(xd)  for  all 
values  of  xd.  Note  that  this  formulation  is  a  generalization  of  the  standard  NLP  problem,  in 
that  it  reduces  to  a  standard  NLP  if  nd  =  0,  in  which  case,  £,  A,  and  u  do  not  change. 

In  solving  an  MVP  problem,  we  note  that  continuous  relaxation  of  the  categorical  variables 
is  not  possible,  which  precludes  the  use  of  branch  and  bound  methods.  Furthermore,  we  assume 
that  the  categorical  variable  space  is  sufficiently  large  such  that  optimizing  over  every  possible 
combination  of  categorical  variable  values  is  intractable.  This  means  that  the  best  we  can 
hope  for  is  a  locally  optimal  solution.  However,  a  notion  of  local  optimality  for  MVP  problems 
is  not  well-defined  because  there  is  no  topology  associated  with  the  nonnumeric  categorical 
variables.  Local  optimality  must  therefore  be  defined  in  terms  of  a  problem-specific  set  of 
discrete  neighbors.  This  is  done  by  defining  a  continuous  set- valued  function  J\f  :  X  — *  2X , 
where  2A  is  the  power  set  of  X  (i.e.,  the  set  of  all  possible  subsets  of  X).  Thus  a  point  y  E  X 
is  a  discrete  neighbor  of  X  if  y  E  M(x).  By  convention,  x  E  A f(x). 

For  a  point  x  E  X  to  be  locally  optimal,  we  mean  that  it  is  locally  optimal  in  the  traditional 
sense  when  the  categorical  variables  are  fixed,  that  f(x)  <  f(y)  for  all  y  E  Af(x),  and  that, 
for  any  y  E  N(x)  if  f(x)  =  f(y),  then  y  is  also  locally  optimal  (in  the  traditional  sense)  when 
categorical  variables  are  fixed.  A  formal  definition  is  given  in  US¬ 
As  an  example,  Kokkolaras,  Audet,  and  Dennis  |22]  studied  the  design  of  a  thermal  insulation 
system  of  fixed  length  consisting  of  a  certain  number  of  insulators  of  various  material  types  and 
thicknesses,  and  a  corresponding  set  of  heat  intercepts  placed  between  each  pair  of  insulators. 
The  heat  intercepts  were  set  at  specified  temperatures,  and  the  objective  was  to  minimize  the 
power  to  keep  one  of  the  ends  at  a  fixed  temperature.  Given  a  specific  design,  a  discrete  neighbor 
was  defined  to  be  any  design  obtained  by  replacing  any  single  insulator  with  one  of  a  different 
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type,  adding  an  insulator  and  corresponding  heat  intercepts  at  any  location,  or  removing  any 
insulator  with  its  adjacent  heat  intercept. 

Pattern  search  is  an  iterative  method  that  generates  a  sequence  of  feasible  points  with  non¬ 
increasing  function  values.  At  each  iteration,  the  objective  function  /  is  evaluated  at  points 
lying  on  a  mesh  in  an  attempt  to  find  a  point  with  a  lower  function  value  than  that  of  the 
incumbent.  The  mesh  is  defined  by  a  set  of  directions  that  form  a  positive  spanning  set  023 
(i.e. ,  a  set  of  directions  such  that  any  vector  in  the  space  can  be  expressed  as  a  nonnegative 
linear  combination  of  these  directions). 

At  each  iteration,  the  mesh  is  the  direct  product  of  the  union  of  a  finite  number  of  lattices 
(2  ^  ^  ^ 
in  Mn  with  the  integer  space  Zn  ,  as  follows.  For  each  combination  i  =  1,2,...,  zmax  of  values 

that  the  categorical  variables  may  take  on,  a  set  of  positive  spanning  directions  Dl  =  GiZi  is 

formed,  where  Gi  £  MnCxnC  is  a  nonsingular  generating  matrix  and  Zj  £  ZnCx\D'\.  The  mesh  is 

then  the  direct  product  of  Xd  with  a  union  of  a  finite  number  of  lattices  in  Xc  centered  at  the 

continuous  part  of  the  current  iterate: 

*max 

Mk  =  |J  Mix  Xd,  Ml  =  [J  {xc  +  XkDlz  :  z  £  zj* M}  C  MnC,  (10) 

i=  1  x£Sk 

where  Ak  >  0  is  the  mesh  size  parameter,  and  Sk  is  the  set  of  all  previously  evaluated  trial 
points  Sk  (with  So  as  the  set  of  initial  points).  The  set  of  discrete  neighbors  is  assumed  to  be 
constructed  such  that  every  neighbor  lies  on  the  current  mesh. 

The  evaluation  of  trial  points  on  the  mesh  is  performed  in  three  steps:  an  optional  SEARCH 
step  and  a  local  POLL  step,  and  an  extended  poll  step.  In  the  search  step,  the  function  is 
evaluated  at  a  finite  number  of  mesh  points  Sk.  The  user  may  specify  virtually  any  finite  strategy 
(including  none)  for  identifying  these  points.  Typically,  this  step  is  more  global  in  nature  and 
can  make  use  of  a  favorite  heuristic,  a  set  of  space- filling  points,  or  specialized  knowledge  the  user 
may  have.  If  /  is  computationally  expensive  to  evaluate,  one  common  approach  is  to  construct 
and  optimize  an  inexpensive  surrogate  function  m  on  the  mesh  at  each  search  step. 

If  the  search  is  unsuccessful  in  finding  an  improved  mesh  point  (i.e.,  one  that  has  a  lower 
function  value  than  the  incumbent),  the  POLL  step  is  performed,  in  which  mesh  points  that  are 
adjacent  to  the  incumbent  are  evaluated,  both  in  the  continuous  sense  (while  holding  the  cate¬ 
gorical  variable  values  constant)  and  in  the  discrete  sense  (the  current  set  of  discrete  neighbors). 

Polling  with  respect  to  continuous  variables  requires  use  of  positive  spanning  sets  in  MnC .  Let 
Dl  denote  the  set  of  poll  directions  corresponding  to  the  i-th  set  of  categorical  variable 
values  for  each  iteration  k.  The  poll  set  Pk{x),  centered  at  a  point  x  £  X ,  is  the  set  of  neighboring 
mesh  points  in  the  directions  Dlk,  while  holding  the  categorical  variables  fixed;  i.e., 

Pk(x)  =  {x}  U  {(xc  +  A kd,  xd)  £  X  :  d  £  D{}  C  Mk.  (11) 

If  polling  with  respect  to  the  continuous  variables  fails  to  find  an  improved  mesh  point  in 
Pk(xk),  polling  is  performed  on  the  current  set  of  discrete  neighbors  A f{xk).  The  objective 
function  /  is  evaluated  at  each  of  the  points  in  Pk{xk)  U  N(xk)  until  a  lower  objective  function 
is  found  or  until  all  these  points  have  been  evaluated. 

If  the  both  the  search  and  poll  steps  fail  to  find  an  improved  mesh  point,  then  the 
extended  poll  step  is  performed.  In  this  step,  additional  polling  is  performed  around  each 


7 


discrete  neighbor  point  y  E  J\f(xk),  whose  objective  function  value  was  only  a  small  amount 
greater  than  that  of  the  incumbent.  This  is  done  with  the  idea  that,  while  y  was  not  an  improved 
mesh  point,  it  is  a  promising  region  of  the  space  in  which  to  search  for  one;  thus  polling  around 
y  may  produce  a  better  objective  function  value  than  the  incumbent.  The  extended  poll  step 
is  performed  whenever  /(y&)  <  f(xk )  +  £&,  where  the  extended  poll  trigger  £*,  satisfies  £&  >  £ 
for  a  fixed  £  >  0.  The  values  for  £&  are  often  set  as  a  percentage  of  the  objective  function  value 
at  the  current  iterate  [4].  Larger  values  of  £  will  result  in  more  extended  polling,  thus  more 
computational  cost  per  iteration,  but  may  result  in  a  better  final  solution. 

For  each  iteration  k  and  for  each  point  yk  E  =  {y  E  Af(xk)  :  f{xk)  <  f(y)  <  f(xk)  +  £&}, 
the  extended  poll  step  begins  by  performing  a  poll  around  the  discrete  neighbor  yk  =  y\  and 
then  continues  polling  around  a  finite  sequence  points  {yk}jLi  satisfying  /(y£)  <  /(y^  x)  for  all 
j  =  2,3, ... ,  Jfc.  The  index  J*,  occurs  when  either  /(y)  <  /(x^)  for  some  y  E  P(y^)  or  until 
/(y)  >  (yj!fe)  for  all  V  £  P{Vkk)-  The  set  of  extended  POLL  points  is  therefore  expressed  as 

/  Jk 

xk(tk)=  U  I  U  pk(yi) 

Vk&A ft  V=1 


If  the  search,  POLL,  or  extended  poll  step  is  successful  in  finding  an  improved  mesh 
point,  that  point  becomes  the  new  incumbent  Xk+i  and  the  mesh  is  retained  or  coarsened.  If  all 
three  steps  fail  to  find  an  improved  mesh  point,  then  the  incumbent  is  retained  (i.e.  Xk+\  =  Xk) 
and  declared  to  be  a  mesh  local  optimizer,  and  the  mesh  is  refined.  Refining  or  coarsening  the 
mesh  is  done  by  updating  the  mesh  size  parameter  according  the  rule, 


^fe+i 


=  r^Afc, 


f  {0, 1, . . .  ,w+}  f{xk+i)<f{xk), 

\  {w;-,  w~  +  1, . . . ,  —  1}  otherwise. 


(13) 


where  r  >  1  is  rational,  and  w~  <  —1  and  w+  >  0  are  integers.  Mesh  coarsening  does  not  affect 
theoretical  convergence  properties.  In  practice,  it  can  slow  convergence,  but  it  can  also  cause 
the  algorithm  to  skip  over  a  local  minimum  and  find  a  better  one  [5J. 

The  full  description  of  the  MVPS  algorithm  is  given  in  Figure  [2j  Convergence  of  the  algo¬ 
rithm  to  suitably  defined  first-order  stationary  points  was  proved  in  |4l  and  HU- 

For  problems  with  linear  constraints,  infeasible  points  are  simply  discarded  without  being 
evaluated  by  the  objective  function.  To  guarantee  that  theoretical  convergence  properties  still 
hold,  the  only  additional  requirement  is  that  the  rule  for  selecting  polling  directions  must  con¬ 
form  to  the  geometry  of  the  nearby  linear  constraint  boundaries  m  H23-  An  algorithm  for 
constructing  conforming  directions  is  given  in  f24]  and  [6]  in  the  nondegenerate  and  degenerate 
cases,  respectively. 


4  Implementation 

Before  presenting  numerical  results,  we  first  describe  how  the  x-ray  tomography  problem  was 
setup  to  be  solved,  specifically  with  respect  to  the  chosen  set  of  discrete  neighbors,  materials 
considered,  data  generation,  and  experimental  conditions. 


Mixed  Variable  Pattern  Search  (MVPS)  Algorithm 

Initialization:  Let  xq  E  X  satisfy  f(x o)  <  oo.  Set  Aq  >  0  and  £  >  0. 


For  k  =  0,1,2,...,  perform  the  following: 


1.  Set  the  extended  poll  trigger  £&  >  £. 

2.  Search  step:  Employ  some  finite  strategy  seeking  an  improved  mesh  point;  i.e., 
xk+i  E  Mk  such  that  f{xk+ 1)  <  f{xk). 

3.  Poll  step:  If  the  search  step  does  not  find  an  improved  mesh  point,  evaluate  / 
at  points  in  Pk(xk )  LiJ\f(xk)  until  an  improved  mesh  point  xk+\  is  found  (or  until  all 
points  have  been  evaluated). 

4.  Extended  Poll  step:  If  the  search  and  poll  steps  does  not  find  an  improved 
mesh  point,  evaluate  /  at  points  in  Xk(£,k)  until  an  improved  mesh  point  xk+\  is 
found  (or  until  all  points  have  been  evaluated). 


Update:  If  search,  poll,  or  extended  pc*^  ^ 
update  Xfc+i,  and  set  Ak+i  >  Ak  according  to  (13), 
otherwise,  set  xk+\  =  xk,  and  set  Ak+i  <  Ak  according  to  (13) 


Figure  2:  MVPS  Algorithm 


4.1  Discrete  Neighbors  and  Materials 

Recall  from  earlier  discussion  that  the  set  of  discrete  neighbors  must  be  defined  by  the  user. 
Thus,  when  a  local  solution  to  an  MVP  problem  is  found,  it  is  always  with  respect  to  this  set. 
Choosing  all  possible  discrete  neighbors  might  result  in  a  better  final  solution,  but  it  would 

^max 

require  \M\k  function  evaluations  at  each  unsuccessful  iteration,  where  \M\  is  the  number 
fc=i 

of  material  types  considered.  This  number  grows  very  large,  even  for  modest  values  of  \M\ 
and  nm ax.  Defining  the  discrete  neighborhood  structure  is  perhaps  the  most  crucial  aspect  of 
our  approach  to  solving  the  x-ray  tomography  problem.  It  must  take  into  account  inherent 
properties  of  the  Abel  transform,  as  well  as  known  properties  of  the  material  types. 

To  construct  this  set,  we  first  introduce  the  idea  of  a  pair  of  materials  being  adjacent  to 
each  other  if  they  share  a  common  property  defined  by  the  user.  In  our  case,  adjacent  materials 
may  be  thought  of  as  materials  that  may  be  confused  with  one  another  in  the  current  object 
configuration.  For  example,  beryllium  and  aluminum  exhibit  similar  properties  with  respect  to 
the  scattering  and  attenuation  of  x-ray  photons  as  they  pass  through  these  materials.  Since 
it  would  be  easy  to  confuse  them  when  analyzing  the  results  of  a  radiograph,  beryllium  and 
aluminum  may  be  considered  adjacent.  However,  since  lead  is  not  easily  confused  with  air,  they 
would  not  be  considered  adjacent  to  each  other. 
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The  adjacency  relationship  among  materials  can  be  represented  as  a  graph,  composed  of 
nodes  that  represent  the  materials  and  edges  that  represent  the  adjacency  relationship  between 
each  pair  of  materials.  The  associated  adjacency  matrix  is  square,  symmetric,  and  binary,  where 
each  row  or  column  represents  a  material  type,  and  a  value  of  1  in  the  ( i,j)th  element  means 
that  material  i  is  adjacent  to  material  j.  To  prevent  redundancy,  a  material  is  not  considered 
adjacent  to  itself,  and  thus  the  adjacency  matrix  will  have  zero  entries  along  the  diagonal. 

Table  [I]  shows  the  complete  list  of  materials  we  considered,  their  abbreviations  (which  are 
used  throughout  this  paper),  and  the  material  class  they  represent.  Seven  of  these  materials 
were  selected  to  be  representative  of  a  larger  group  of  materials  that  all  possess  similar  physical 
attributes  and  interact  with  x-ray  photons  in  a  similar  manner.  Copper  (Cu)  and  Iron  (Fe)  were 
deliberately  chosen  to  represent  very  similar  material  groups  in  order  to  test  the  algorithm’s 
ability  to  distinguish  between  two  different  materials  with  very  similar  physical  attributes. 

The  adjacency  matrix  we  used  is  given  in  Table  [2]  It  represents  a  connected  graph,  so  that 
any  two  materials  can  be  reached  from  each  other  through  a  finite  set  of  adjacency  relationships. 
This  is  a  key  property  in  the  optimization  process  [TJ ,  as  it  (in  theory)  allows  all  possible  material 
configurations  to  be  considered. 


Table  1:  Materials  Considered. 


Material 

Abbreviation 

Representative  class 

Air 

Air 

Voids  and  unpressurized  gases 

Polyethylene 

Peth 

Light  plastics  and  organic  materials 

Beryllium 

Be 

Beryllium 

Teflon 

Teflon 

Dense  plastics 

Aluminum 

A1 

Aluminum 

Stainless  Steel 

Fe 

Iron-like  medium  density  metals 

Copper 

Cu 

Other  medium  density  metals 

Lead 

Pb 

Heavy  metals 

Uranium 

U 

Very  heavy  metals 

Given  a  specific  configuration,  we  included  several  different  types  of  discrete  neighbors, 
obtained  by  doing  one  of  the  following: 

1.  Swap  any  layer  with  a  material  of  a  different  type, 

2.  Delete  any  single  material  layer, 

3.  Insert  a  single  material  layer  between  any  two  existing  layers  or  at  the  object  center, 

4.  Split  any  layer  into  a  pair  of  adjacent  layers  of  different  types, 

5.  Merge  any  two  layers  together, 

6.  Remove  simultaneously  all  layers  with  less  than  a  specified  minimum  thickness, 

7.  Combine  all  physically  adjacent  layers  of  the  same  material  type. 

We  limit  this  set  of  neighbors  using  the  adjacency  matrix,  along  with  some  other  reasonable 
restrictions.  In  describing  each  in  greater  detail,  we  index  the  layers  i  =  1,2,  ...,n,  ordered 
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Table  2:  Adjacency  Matrix  for  Materials. 


Air 

Peth 

Be 

Teflon 

A1 

Fe 

Cu 

Pb 

u 

Air 

0 

1 

1 

1 

1 

0 

0 

0 

0 

Peth 

1 

0 

1 

1 

1 

0 

0 

0 

0 

Be 

1 

1 

0 

1 

1 

0 

0 

0 

0 

Teflon 

1 

1 

1 

0 

1 

1 

1 

0 

0 

A1 

1 

1 

1 

1 

0 

1 

1 

0 

0 

Fe 

0 

0 

0 

1 

1 

0 

1 

1 

1 

Cu 

0 

0 

0 

1 

1 

1 

0 

1 

1 

Pb 

0 

0 

0 

0 

0 

1 

1 

0 

1 

U 

0 

0 

0 

0 

0 

1 

1 

1 

0 

from  interior  to  exterior.  The  transmission  layer  n  +  1  has  fixed  material  type  mn. |_i  and  edge 
location  xn+i  =  L. 

A  swap  of  materials  in  a  layer  involves  replacing  the  material  of  one  layer  with  a  different 
type  from  the  material  library  M  while  holding  all  layer  thicknesses  fixed.  The  new  material 
must  be  adjacent  to  the  old  material  (as  defined  by  the  adjacency  matrix),  and  mn  /  mn+ 1 
must  always  hold. 

Deletion  of  any  layer  i  (with  its  material  mi  and  edge  location  Xi)  may  be  performed  as  long 
as  the  number  of  layers  n  does  not  fall  below  nmin.  When  layer  i  is  deleted,  either  layer  i  —  1  or 
i  +  1  can  be  chosen  to  take  on  the  thickness  of  the  deleted  layer  and  to  become  layer  i.  When 
i  =  n  >  nm in  +  2  (i.e. ,  the  outermost  layer)  is  considered  for  deletion,  and  layer  n  —  1  is  of  the 
same  material  type  as  the  transmission  layer  (i.e.,  mn_ i  =  mn+i),  then  layers  n  and  n  —  1  are 
both  deleted  at  once. 

Insertion  of  an  additional  layer  between  any  two  layers  i  and  i  +  1  may  be  done  with  any 
material  that  is  adjacent  (as  defined  by  the  adjacency  matrix)  to  either  m*  or  mi+ 1,  provided 
that  n  does  not  exceed  nmax  after  insertion  and  provided  layers  i  and  i  + 1  are  sufficiently  thick. 
In  this  case,  the  new  edge  location  is  set  between  Xi  and  Xj+i:  either  xt  is  reduced  enough  to 
allow  the  new  layer’s  edge  location  to  take  the  old  value  of  xt  (the  inner  layer’s  thickness  is 
reduced),  or  else  the  new  edge  location  is  set  to  be  far  enough  between  xt  and  Xj+i  so  that  Xi 
can  remain  constant  (the  outer  layer’s  thickness  is  reduced). 

Splitting  a  layer  not  only  divides  the  layer  exactly  in  half,  but  also  changes  each  material 
type  to  an  adjacent  one  (actually,  the  next  adjacent  one  in  the  list  of  adjacent  materials).  For 
example,  A1  would  split  into  Teflon  and  Fe,  in  either  order. 

Merging  two  layers  occurs  when  two  layers  are  combined  and  the  material  type  of  the  new 
layer  is  chosen  to  be  the  weighted  average  of  the  material  types,  relative  to  the  indexing  of  the 
adjacency  matrix.  For  example,  Teflon  and  Peth  would  be  merged  into  Be. 

For  neighbors  having  all  thin  layers  removed  simultaneously,  either  the  interior  or  exterior 
physically  adjacent  layer  is  expanded  to  compensate. 

Finally,  since  the  order  in  which  points  are  evaluated  can  have  a  significant  effect  on  the 
solution,  two  ordering  strategies  were  employed.  First,  neighbors  that  result  in  fewer  layers 
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(deletion,  merging,  etc.)  were  evaluated  first,  while  those  that  increase  the  number  of  layers 
(insertion,  splitting,  etc.)  were  evaluated  last.  Second,  because  the  outside  layers  affect  the  x- 
rays  more  than  the  inside  layers  do,  neighbors  with  changes  to  outside  layers  are  always  evaluated 
before  those  with  changes  to  internal  layers. 

4.2  Experimental  Conditions  and  Data  Generation 

The  x-ray  beam  was  modeled  as  produced  by  a  2.4  MeV  Cygnus  source.  We  approximated  this 
continuous  source  as  23  discrete  monoenergetic  sources  of  weighted  intensities.  Energies  varied 
from  100  keV  to  2.3  MeV  in  100  keV  intervals.  The  beam  geometry  for  all  experiments  was 
plane  parallel  with  uniform  cross-sectional  intensity. 

The  object  symmetry  axis  and  detector  plane  were  oriented  perpendicular  to  the  incident 
beam  propagation  direction.  In  this  paper  we  do  not  model  scattering,  so  the  distances  between 
source,  object,  and  detector  have  no  effect  on  the  results.  The  detector  is  given  a  resolution  of 
w  =  200  jan.  and  has  100%  energy  collection  efficiency. 

Simulated  radiographs  d  are  produced  by  constructing  a  model  of  the  forward  projection 
operator  P.  We  used  two  distinct  methods,  both  which  utilize  the  Abel  transform.  A  linear 
polychromatic  (LP)  approximation  is  the  superposition  of  linear  attenuation  radiographs  at 
each  of  the  23  distinct  energy  levels  of  the  Cygnus  source.  A  linear  monochromatic  (LM) 
approximation  uses  a  single  attenuation  scalar  for  each  material  that  is  the  source-intensity 
weighted  average  of  attenuation  coefficients  for  each  material.  Both  methods  partially  account 
for  the  energy  dependence  of  the  attenuation  coefficients  but  in  different  ways.  Simulated  Poisson 
noise  was  added  to  the  radiograph  data  assuming  30,000  photons  per  pixel  (unattenuated).  This 
provides  a  realistic  but  technically  inaccurate  noise  profile. 

5  Numerical  Testing 

To  demonstrate  the  effectiveness  of  the  MVPS  algorithm  in  quantitatively  reconstructing  objects, 
we  applied  the  NOMADrn  MATLAB®  software  package  |3'  to  three  test  sets  that  differ  in 
their  construction,  with  each  set  consisting  of  6-10  test  object  configurations.  Following  a  brief 
description  of  test  conditions  and  parameter  settings,  we  describe  each  test  set  in  turn  and 
present  numerical  results. 

MVPS  parameters  were  set  as  follows.  The  extended  poll  trigger  was  set  at  0.99.  At  this 
value,  extended  polling  occurred  whenever  a  neighboring  object’s  configuration  resulted  in  an 
objective  function  value  that  was  less  than  99%  higher  than  that  of  the  current  best  configuration. 
The  initial  mesh  size  was  set  at  0.1  cm  and  doubled  (r  =  2,  Wk  =  1)  whenever  an  improved  mesh 
point  was  found,  provided  the  mesh  size  did  not  exceed  0.1  cm.  If  not  found,  the  mesh  size  was 
divided  in  half  (u>k  =  —  1).  Termination  occurred  the  first  time  the  mesh  size  fell  below  10-3. 
No  search  step  was  used  and  the  poll  step  used  the  standard  coordinate  directions  and  their 
negatives. 

To  keep  computational  time  reasonable,  the  minimum  and  maximum  number  of  allowable 
material  layers  were  set  to  nmin  =  1  and  ?rmax  =  6,  respectively,  and  the  minimum  layer  thickness 
was  set  at  5  =  0.08  cm.  We  assume  that  we  can  independently  obtain  a  good  estimate  of  the 
object  radius  r,  and  we  set  the  initial  guess  to  be  a  cylinder  of  Al,  but  represented  as  two  adjacent 
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layers  of  equal  thickness;  i.e. ,  x  =  [r/2,r]  and  m  =  { Al,  Al}.  The  representation  of  two  layers  of 
the  same  material  is  done  to  help  the  algorithm  avoid  early  termination  due  to  limitations  of  the 
discrete  neighbor  set.  Aluminum  was  chosen  as  the  initial  material  guess  because  it  represents 
the  “center”  of  the  material  library. 

5.1  Test  Scenarios 

Scenario  1.  The  first  test  set  consists  of  six  objects  with  3  cm  radii,  each  having  three  layers 
whose  material  types  are  the  six  possible  orderings  of  Fe,  Be,  and  Peth,  and  whose  thickness 
are  all  equal  to  1  cm  (See  Table  [3]) .  Radiographs  were  computed  using  the  LM  method.  The 
material  library  used  in  the  MVP  algorithm  included  all  the  materials  from  Table  [l]  except  Cu. 
The  adjacency  matrix  was  that  of  Table  [2]  with  the  row  and  column  for  Cu  removed.  The  cost 
function  computes  a  relative  merit  value  using  the  same  LM  approximation  used  in  radiograph 
construction. 

This  experiment  is  ideal  for  initial  testing  because  the  material  types  are  distinct  in  their 
attenuation  characteristics  and  because  the  cost  function  and  simulation  utilize  the  same  forward 
measurement  operator  P,  allowing  for  the  possibility  of  an  exact  solution,  modulo  the  effects  of 
data  noise. 


Table  3:  Test  Set  1:  Object  Configurations 


Object 

Materials 

Edge  Locations  (cm) 

la 

[Fe,  Be,  Peth] 

[1.0,  2.0,  3.0] 

lb 

[Fe,  Peth,  Be] 

[1.0,  2.0,  3.0] 

lc 

[Be,  Fe,  Peth] 

[1.0,  2.0,  3.0] 

Id 

[Be,  Peth,  Fe] 

[1.0,  2.0,  3.0] 

le 

[Peth,  Fe,  Be] 

[1.0,  2.0,  3.0] 

If 

[Peth,  Be,  Fe] 

[1.0,  2.0,  3.0] 

Scenario  2.  The  second  test  set  consists  of  10  randomly  generated  configurations  of  1-4 
layers  with  an  enforced  minimum  layer  thickness  of  0.1  cm.  The  possible  object  materials  were 
limited  to  a  subset  of  the  full  material  library:  Air,  Peth,  Be,  Teflon,  Al,  and  Fe.  These  objects 
are  described  in  Table  [4]  The  radiographs  were  computed  using  the  LP  method.  The  material 
library,  adjacency  matrix,  and  cost  function  used  by  the  MVP  algorithm  were  the  same  as 
Scenario  1. 

This  experiment  was  used  to  test  the  ability  of  the  algorithm  to  identify  objects  of  greater 
variety  of  materials  and  layer  thicknesses.  It  was  also  used  to  examine  the  effects  of  using  a  fully 
linearized  measurement  operator  (LM)  in  the  cost  function  for  examining  radiographs  produced 
by  considering  approximate  polychromatic  effects  (LP). 

Scenario  3.  The  final  scenario  is  identical  to  Scenario  2  except  for  the  cost  function,  which 
utilizes  the  LP  method.  So  as  in  Scenario  1,  the  cost  function  and  simulation  utilize  the  same 
measurement  operator  allowing  for  the  possibility  of  an  exact  solution,  modulo  the  effects  of 
data  noise. 
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Table  4:  Test  Set  2:  Object  Configurations. 


Object 

Materials 

Edge  Locations  (cm) 

2a 

[Fe,  Teflon,  Fe] 

[1.3973,  1.7028,  2.7225] 

2b 

[Be] 

[2.4927] 

2c 

[Be,  Fe] 

[1.7136,  2.4441] 

2d 

[Air,  Al] 

[3.6109,  3.8152] 

2e 

[Be,  Air,  Teflon,  Be] 

[0.2379,  0.8433,  1.9587,  3.4136] 

2f 

[Al,  Teflon,  Be] 

[0.2151,  2.6260,  3.7330] 

2g 

[Fe,  Be] 

[1.4035,  2.8271] 

2h 

[Peth,  Al,  Air,  Fe] 

[1.7489,  2.4271,  3.0819,  3.6769] 

2i 

[Air,  Al,  Teflon,  Al] 

[0.8161,  1.8739,  2.0518,  3.7236] 

2j 

[Be,  Peth,  Fe] 

[1.3628,  1.9025,  3.6278] 

5.2  Results 

For  each  object,  a  run  was  made  in  which  all  iterates  were  recorded,  and  the  best  point  for  each 
setting  of  the  categorical  variable  values  was  extracted,  with  the  restriction  that  any  solution 
with  any  layer  width  less  than  0.1  cm  was  discarded.  This  process  led  to  a  small  set  of  meaningful 
candidate  solutions.  Since  we  actually  know  the  true  solution  for  each  object,  we  performed  an 
additional  run  from  a  different  initial  point  whenever  we  failed  to  find  a  point  whose  function 
value  came  within  5%  of  that  of  the  true  solution.  The  second  initial  point  is  the  same  as  the 
first,  except  the  inner  material  is  changed  to  Air;  i.e.,  x  =  [r/2,r],  m  =  {Air,  Al} 

In  describing  the  results  that  follow,  we  denote  by  nSol  the  total  number  of  solutions  found 
whose  function  value  was  within  5%  of  that  of  the  true  solution,  and  we  denote  by  rank  the 
rank  of  the  “true”  solution  in  the  list  of  best  points  found.  For  each  object  we  attach  a  rating 
of  H  (Hit),  I  (Inconclusive),  or  M  (Miss),  which  has  the  following  interpretation: 

1.  H  ( nSol  >  0  and  rank  <  nSol ):  This  means  that  there  is  a  set  of  solutions  and  the  true 
solution  is  among  them.  Ideally,  we  want  nSol  =  rank  =  1,  but  we  recognize  that  a  small 
set  of  ranked  solutions  is  often  a  more  meaningful  result. 

2.  I  (nSol  =  0):  This  is  the  restart  indicator;  none  of  the  solutions  seem  to  adequately 
describe  the  data,  from  which  we  conclude  that  we  are  stuck  at  a  local  minimizer. 

3.  M  ( nSol  >  0  and  the  true  solution  is  unranked  or  not  present):  This  is  the  worst  situation 
because  there  appear  to  be  good  solutions,  so  that  a  restart  is  not  indicated,  but  the  true 
solution  is  not  in  the  solution  set. 

Our  motivation  for  reporting  results  in  this  manner  is  as  follows.  Recall  that  the  objective 
function  value  of  a  solution  is  a  measure  of  how  accurately  the  model  and  solution  together 
describe  the  data  (see  <§)■  An  objective  function  value  of  zero  would  indicate  a  perfect  match, 
but  even  the  true  solution  does  not  achieve  this  because  of  the  noise  in  the  data.  Because  of  both 
model  deficiencies  and  data  noise,  it  makes  sense  to  consider  candidate  solutions  with  roughly 
the  same  objective  function  value  as  the  true  solution;  i.e.,  we  do  not  want  to  arbitrarily  dismiss 
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good  local  solutions.  The  choice  of  5%  in  the  discussion  above  was  somewhat  arbitrary,  but 
the  results  are  robust  to  the  particular  value.  When  the  model  is  reasonably  accurate,  a  true 
solution  in  categorical  variables  should  attain  an  objective  function  value  well  within  the  5% 
condition. 

Tables  id  summarize  the  results  for  the  three  Test  Sets,  along  with  reruns  for  the  last  two 
scenarios.  Solutions  were  obtained  in  a  matter  of  minutes  on  a  single  processor  PC,  which  was 
consistent  with  our  goal  of  achieving  solutions  quickly,  especially  for  time-sensitive  applications. 
The  median  and  mean  solution  times  were  10  and  12  minutes,  respectively,  and  the  maximum 
solution  time  was  33  minutes  for  the  initial  run  of  Problem  2h. 


Table  5:  Scenario  1  results  show  an  excellent  6-0-0  (Hit-Inconclusive-Miss)  score. 


Object 

nSol 

rank 

Rating 

la 

1 

1 

H 

ib 

1 

1 

H 

lc 

5 

1 

H 

id 

4 

2 

H 

le 

6 

2 

H 

if 

1 

1 

H 

Table  6:  Scenario  2  results  score  a  4-5-1  (Hit-Inconclusive- Miss),  but  improve  to  6-3-1  when  we 
include  results  obtained  from  the  second  starting  point. 


Object 

Initial  run 

Restart 

nSol 

rank 

Rating 

nSol 

rank 

Rating 

2a 

5 

1 

H 

5 

1 

H 

2b 

1 

1 

H 

1 

1 

H 

2c 

2 

- 

M 

2 

- 

M 

2d 

0 

2 

I 

1 

1 

H 

2e 

5 

1 

H 

5 

1 

H 

2f 

2 

1 

H 

2 

1 

H 

2g 

0 

1 

I 

0 

1 

I 

2h 

0 

1 

I 

1 

1 

H 

2i 

0 

- 

I 

0 

- 

I 

2j 

0 

- 

I 

0 

- 

I 

Table  [5] shows  excellent  performance  for  the  scenario  1  set  of  test  problems.  In  all  cases,  the 
initial  point  led  to  a  set  of  solutions  that  contained  the  true  solution.  The  largest  solution  set 
contained  six  candidates.  In  no  case  was  the  true  solution  ranked  worse  than  second.  Figure  [3] 
illustrates  the  rating  of  H  (hit)  for  object  Id.  The  true  object  (a)  and  the  first  three  candidate 
solutions  (b-d)  are  shown  along  with  the  corresponding  objective  function  values  /.  All  of 
the  objects  in  the  reconstruction  set  have  similar  objective  function  values.  The  optimization 
method  returns  very  good  solutions,  but  the  objective  function  itself  has  difficulty  distinguishing 
objects  differing  in  categorical  variables.  All  of  the  reconstructions  shown  provide  equally  likely 
explanations  of  the  data. 


15 


Table  7:  Scenario  3  results  score  a  6-4-0  (Hit-Inconclusive- Miss),  but  improve  to  9-1-0  when  we 
include  results  obtained  from  the  second  starting  point. 


Object 

Initial  run 

Restart 

nSol 

rank 

Rating 

nSol 

rank 

Rating 

3a 

4 

2 

H 

4 

2 

H 

3b 

1 

1 

H 

1 

1 

H 

3c 

2 

2 

H 

2 

2 

H 

3d 

0 

2 

I 

1 

1 

H 

3e 

0 

- 

I 

7 

4 

H 

3f 

2 

1 

H 

2 

1 

H 

3g 

1 

1 

H 

1 

1 

H 

3h 

1 

1 

H 

1 

1 

H 

3i 

0 

- 

I 

0 

- 

I 

3j 

0 

1 

I 

1 

1 

H 

In  Table  [6]  we  note  decent  performance  for  the  scenario  2  set  of  test  problems.  Figure  [4] 
illustrates  the  rating  of  M  (miss)  for  object  2c.  The  true  object  (a)  and  the  first  three  candidate 
solutions  (b-d)  are  shown  along  with  the  corresponding  objective  function  values  /.  The  first 
two  reconstructions  attain  objective  function  values  with  the  5%  criterion  and  the  third  is 
similar,  yet  the  true  solution  (in  categorical  variables)  is  not  present.  For  the  entire  test  set, 
we  obtain  an  initial  H-I-M  result  of  4-5-1,  which  improves  to  6-3-1  with  the  addition  of  the 
five  restart  results.  However,  some  problems  could  not  be  solved  adequately.  We  attribute  this 
mediocre  performance  to  model  inconsistency.  The  cost  function  uses  a  monochromatic  beam 
approximation  while  the  simulated  radiographs  were  computed  using  a  linear  polychromatic 
approximation.  In  spite  of  this,  we  failed  to  generate  a  solution  ( nSol  >  0)  for  only  one  test 
problem. 

In  Table  [7]  we  note  excellent  performance  for  the  scenario  3  set  of  test  problems.  We  are 
able  to  solve  9  of  10  problems  (using  both  initial  points),  and  the  remaining  problem  (3i)  still 
indicates  a  restart.  No  missed  solutions  are  indicated.  Figure  [5]  illustrates  the  rating  of  I 
(inconclusive)  for  object  3e  before  restart.  The  true  object  (a)  and  the  first  three  candidate 
solutions  (b-d)  are  shown  along  with  the  corresponding  objective  function  values  /.  No  solution 
provides  a  reasonable  objective  function  value.  This  indicates  that  these  solutions  reside  in  very 
nonoptimal  local  minima. 

6  Conclusions 

Quantitative  reconstruction  of  cylindrically  symmetrical  objects  from  x-ray  data  normally  re¬ 
quires  an  excessive  amount  of  time  to  accomplish,  and  fast  methods,  such  as  regularized  inver¬ 
sions  fail  to  produce  a  quantitative  description.  This  work  represents  the  first  real  success  in 
this  area,  in  which  the  result  is  a  quantitative  description  of  the  object  computed  in  a  reasonable 
amount  of  time.  MVPS  is  able  to  rapidly  reconstruct  an  unknown  object’s  composition  from 
x-ray  radiographs  because  of  its  ability  to  handle  mixed  variable  problems  and  incorporate  prior 
knowledge  about  the  object. 
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Figure  3:  H  (hit)  rating  example  for  test  Id  showing  the  true  object  (a),  the  best  three  re¬ 
constructions  (b-d),  and  corresponding  objective  function  values  /.  Radial  positions  of  layer 
boundaries  are  given  in  centimeters  and  material  identifications  are  provided  within  or  near 
each  layer. 


Figure  4:  M  (miss)  rating  example  for  test  2c  showing  the  true  object  (a),  the  best  three 
reconstructions  (b-d),  and  corresponding  objective  function  values  /.  Radial  positions  of  layer 
boundaries  are  given  in  centimeters  and  material  identifications  are  provided  within  or  near  each 
layer. 
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Figure  5:  I  (inconclusive)  rating  example  for  test  3e  showing  the  true  object  (a),  the  best  three 
reconstructions  (b-d),  and  corresponding  objective  function  values  /.  Radial  positions  of  layer 
boundaries  are  given  in  centimeters  and  material  identifications  are  provided  within  or  near  each 
layer. 


Clearly,  our  approach  is  not  foolproof.  Some  objects  are  very  challenging  to  reconstruct 
because  trial  points  are  sensitive  to  either  small  changes  in  edge  locations  or  the  initial  point. 
Including  more  discrete  neighbors  does,  in  fact,  tend  to  improve  solutions,  but  it  also  drives 
up  the  computational  cost.  Developing  strategies  to  improve  this  process  remains  an  open  area 
of  research.  The  success  of  our  approach  relies  on  the  existence  of  a  material  library  that  is 
reasonably  short  and  which  contains  all  materials  likely  to  be  found  in  the  object  of  interest. 

This  approach  to  x-ray  tomography  is  not  limited  to  cylindrically  symmetric  objects.  Our 
formulation  is,  in  a  generalized  sense,  an  optimization  problem  in  a  space  of  categorical  variables 
(materials  to  consider)  and  continuous  variables  (object  geometric  description).  The  application 
potential  for  parametrized  object  identification  from  radiographs  is  clear.  This  approach  is  also 
not  limited  to  particular  assumptions  on  problem  linearity.  The  measurement  operator,  P, 
can  be  any  appropriate  forward  projection  and  can  include  such  diverse  information  as  source 
spectrum,  beam  geometry,  nonlinear  photon/matter  interaction,  detector  response,  and  data 
postprocessing. 
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